getwd()
setwd(getwd())
#setwd("D:\thasegawa\ExtremeEvent\R")
library(reshape2)
library(ggplot2)
library(plyr)
library(grid)
library(scales)

regionlabel <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="USA"] <- "United States"
  levels(y$Region)[levels(y$Region)=="XE25"] <- "EU25"
  levels(y$Region)[levels(y$Region)=="XER"] <- "Rest of Europe"
  levels(y$Region)[levels(y$Region)=="TUR"] <- "Turkey"
  levels(y$Region)[levels(y$Region)=="XOC"] <- "Oceania"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="JPN"] <- "Japan"
  levels(y$Region)[levels(y$Region)=="XSE"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="XSA"] <- "Rest of South Asia"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="XLM"] <- "Rest of South America"
  levels(y$Region)[levels(y$Region)=="CIS"] <- "FSU"
  levels(y$Region)[levels(y$Region)=="XME"] <- "Middle East"
  levels(y$Region)[levels(y$Region)=="XNF"] <- "North Africa"
  levels(y$Region)[levels(y$Region)=="XAF"] <- "Sub-Sa Africa"
z <- y
return(z)
}

croplabel <- function(y){
  levels(y$Crop)[levels(y$Crop)=="PDR"] <- "rice"
  levels(y$Crop)[levels(y$Crop)=="WHT"] <- "wheat"
  levels(y$Crop)[levels(y$Crop)=="GRO"] <- "other grains"
  levels(y$Crop)[levels(y$Crop)=="OSD"] <- "oil crops"
  z <- y
  return(z)
}
fsize <- 10
MyThemeLine <- theme_grey() +
  theme(
    panel.border=element_rect(fill=NA),
    title=element_text(size=fsize,face="bold"),
    axis.title=element_text(size=fsize,face="bold"),
    axis.text = element_text(hjust=1,size = fsize, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.background=element_blank(),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=fsize, colour = "black", angle = 0,face="bold"),
    legend.title = element_text(size=fsize),
    legend.text = element_text(size=fsize)
  )

#Defining function
#Making two figures
twofigure <- function(p1,p2,p3,p4,x,l,wd,ht,rs,x2){
  grid.newpage()
  png(filename=file_name_out, width=wd, height=ht,res=rs)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=3,layout.pos.col=1))
  print(p4, vp = viewport(layout.pos.row=3,layout.pos.col=2))
  #  print(title=x2)
  popViewport()
  dev.off()  
  return()
}

colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","Crop","Value")
pdata <- read.table('../data/dataall.txt')
names(pdata) <- colnames

pdata$RCP <- factor(pdata$RCP, levels=c("RCP85","RCP26"))
pdata$Region <- factor(pdata$Region, levels=c("USA","CAN","JPN","TUR","XER","XE25","XOC","CIS","XNF","CHN","BRA","IND","XLM","XSE","XME","XSA","XAF","WLD"))

#---- making directory
dirname0 <- paste('../output/png/', sep = '')
dir.create(dirname0)

var.in <- c("PCAL","PPOP_UNSH","YEXO","XPRP")
var.in <- c("PCAL")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  dirname <- paste(dirname0, varname, sep = '')
	dir.create(dirname)

}

#-----PPOP_UNSH & PCAL: Histogram for a single region ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
var.in <- c("PCAL","PPOP_UNSH")
regionname <- c("WLD")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
if (varname=="PCAL"){
  labelname <- c("calorie intake [kcal/person/day]")
} else if (varname=="PPOP_UNSH"){
  labelname <- c("population at risk of hunger [billion]")
} else if (varname=="YIELD"){
  labelname <- c("yield [tonne/ha]")
} else {
  labelname <- c("aaaaaaaa")
}

#filenamehist <- paste('../output/png/', varname ,'/',regionname,'_',yearname,'_limited.png', sep = '')
filenamehist <- paste('../output/png/', varname ,'/',regionname,'_',yearname,'.png', sep = '')

xhist <- pdata[pdata$VAR==varname & pdata$Year==yearname & pdata$Region==regionname  & pdata$SSP=="SSP2" & pdata$Crop=="dum", c(-1,-2,-3,-4,-10)]
nocc1 <- pdata[pdata$VAR==varname & pdata$Year==yearname & pdata$Region==regionname & pdata$SSP=="SSP2" & pdata$GCM=="ge2m" & pdata$RCP=="RCP85" & pdata$SYR==yearname & pdata$CO2fert=="noco2" & pdata$N=="N1"& pdata$Crop=="dum", c(-1,-2,-3,-4,-5,-6,-7,-8,-9,-10,-11)]

names(xhist) <- c("GCM","RCP","CO2fert","SYR","N","Value","NoCC")	#Rename the column

minv <- min(xhist[6],xhist[7])

if (varname=="PPOP_UNSH"){
#  maxv <- c(0.15)
  maxv <- max(xhist[6],xhist[7])
    bw <- maxv/100
  } else {
  maxv <- max(xhist[6],xhist[7])
  bw <- diff(range(xhist$Value))/100
}

mainlabel <- c(regionname,yearname)
interact<-interaction(xhist$RCP, xhist$CO2fert, sep=" : ")

mu <- ddply(xhist, "interact", summarise, grp.mean=mean(Value))
v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))

ghist <- ggplot(data=xhist, aes(x=Value,fill=interact)) +
  geom_histogram(alpha = 0.3, position = "identity", aes(xhist=..density..), binwidth=bw, color="black") +
#        geom_density(alpha = 0.3, adjust = 0.2, position = "identity", stat = "density") +
#        geom_density(alpha = 0, adjust = 0.2, position = "identity", stat = "density") +
  xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
  ggtitle(mainlabel) +
  ylab("density") + xlab(labelname) +
  MyThemeLine +
	geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
#  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
#  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
# geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1.2) 
#  geom_text(aes(x=nocc1,y=0.023,label=nocc1),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001)

#  geom_area(mapping = aex(y = ifelse(y<0.2)))

plot(ghist)

ggsave(ghist,filename=filenamehist, dpi=250,width=8, height=6)

}

#----End --- PPOP_UNSH & PCAL: Histogram for a single region----


#----- PPOP_UNSH & PCAL: Boxplot across regions ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
var.in <- c("PCAL","PPOP_UNSH")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  
if (varname=="PCAL"){
  labelname <- c("calorie intake [kcal/person/day]")
} else if (varname=="PPOP_UNSH"){
  labelname <- c("population at risk of hunger [billion]")
} else if (varname=="YIELD"){
  labelname <- c("yield [tonne/ha]")
} else {
  labelname <- c("aaaaaaaa")
}

#  filenamebox <- paste('../output/png/', varname ,'/',yearname,'_limited.png', sep = '')
  filenamebox <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')

#colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")

#	y <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Crop=="dum" & pdata$Region!="WLD" & 
#               pdata$Region!="CHN" & pdata$Region!="BRA" & pdata$Region!="CIS" & pdata$Region!="XNF", c(-1,-2,-4,-10)]

  y <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Crop=="dum" & pdata$Region!="WLD", c(-1,-2,-4,-10)]


	names(y) <- c("Region","GCM","RCP","CO2fert","SYR","N","Value","NoCC")	#Rename the column
  interact<-interaction(y$RCP, y$CO2fert, sep=" : ")

	minv <- min(y[7],y[8])
	maxv <- max(y[7],y[8])

#	bw <- diff(range(y$Value))/50
#	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
#	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))


  y0 <- ddply(y,c("Region","interact"),transform,v0=min(Value))
  y1 <- ddply(y0,c("Region","interact"),transform,v25=quantile(Value,0.25))
  y2 <- ddply(y1,c("Region","interact"),transform,v50=median(Value))
  y3 <- ddply(y2,c("Region","interact"),transform,v75=quantile(Value,0.75))
  y4 <- ddply(y3,c("Region","interact"),transform,v100=max(Value))
	z <- regionlabel(y4)

gbox <- ggplot() + 
    geom_boxplot(data=z, aes(x=Region,y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
  	ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
#	ggtitle(mainlabel)
	ylab(labelname) + xlab("") +	MyThemeLine +
  geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")

  gbox1 <- gbox + coord_flip()
#  g2 <- g1 + facet_wrap(~RCP,ncol=2)

plot(gbox1)

ggsave(gbox1,filename=filenamebox, dpi=250,width=10, height=6)

}

#----End --- PPOP_UNSH & PCAL: Boxplot across regions ----

#-----  YIELD(YEXO) & price(XPRP): Boxplot across regions ----


var.in <- c("XPRP")
var.in <- c("YEXO")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
  } else if (varname=="XPRP"){
    labelname <- c("production price [2005=1]")
  } else {
    labelname <- c("aaaaaaaa")
  }
  
#  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'_limited.png', sep = '')
  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')
  
  #colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")
  
  # Limit the num of regions
  #pyld <- pdata[pdata$VAR=="YEXO" & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD" & 
  #               pdata$Region!="XOC" & pdata$Region!="XNF" & pdata$Region!="XER" & pdata$Region!="XE25" & 
  #               pdata$Region!="USA" & pdata$Region!="TUR" & pdata$Region!="JPN" & pdata$Region!="CIS" & 
  #               pdata$Region!="CAN"& pdata$Region!="BRA" & pdata$Region!="CHN", c(-1,-2,-4)]
  
  # Limit the num of regions
  pyld <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD", c(-1,-2,-4)]
  
  
  names(pyld) <- c("Region","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")	#Rename the column
  
  minv <- min(pyld[8],pyld[9])
  maxv <- max(pyld[8],pyld[9])
#  maxv <- c(20)
  
  #	bw <- diff(range(y$Value))/50
  #	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
  #	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))
  
  interact<-interaction(pyld$RCP, pyld$CO2fert, sep=" : ")
  
  pyld0 <- ddply(pyld,c("Region","interact","Crop"),transform,v0=min(Value))
  pyld1 <- ddply(pyld0,c("Region","interact","Crop"),transform,v25=quantile(Value,0.25))
  pyld2 <- ddply(pyld1,c("Region","interact","Crop"),transform,v50=median(Value))
  pyld3 <- ddply(pyld2,c("Region","interact","Crop"),transform,v75=quantile(Value,0.75))
  pyld4 <- ddply(pyld3,c("Region","interact","Crop"),transform,v100=max(Value))
  pyld5 <- regionlabel(pyld4)
  pyldz <- croplabel(pyld5)
  
  
  gyldbox <- ggplot() + 
    geom_boxplot(data=pyldz, aes(x=Crop, y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
    ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) + ggtitle("") +
    ylab(labelname) + xlab("") +	MyThemeLine +
    geom_point(data=pyldz, aes(x=Crop, y=NoCC),size=6,shape=73,colour="black")
  
  gyldbox1 <- gyldbox + coord_flip()
  gyldbox2 <- gyldbox1 + facet_wrap(~Region, ncol=4)
  
  plot(gyldbox2)
  
  ggsave(gyldbox2,filename=filenamebox2, dpi=250,width=7.5, height=7.5)
  
}



#----- END --- YIELD : Boxplot across regions ----


#----- PPOP_UNSH & PCAL: Bars across regions ----

var.in <- c("PPOP_UNSH")
var.in <- c("PCAL")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
if (varname=="PCAL"){
  labelname <- c("calorie intake [kcal/person/day]")
} else if (varname=="PPOP_UNSH"){
  labelname <- c("population at risk of hunger [billion]")
} else if (varname=="YIELD"){
  labelname <- c("yield [tonne/ha]")
} else {
  labelname <- c("aaaaaaaa")
}

name <- paste('../output/png/', varname ,'/',yearname,'_stack_limit.png', sep = '')
#name <- paste('../output/png/', varname ,'/',yearname,'_stack.png', sep = '')

y <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Crop=="dum", c(-1,-2,-4,-10)]

names(y) <- c("Region","GCM","RCP","CO2fert","SYR","N","Value","NoCC")	#Rename the column
interact<-interaction(y$RCP, y$CO2fert, sep=" : ")

y5 <- ddply(y,c("Region","interact"),transform,v1=quantile(Value,0.01))
y6 <- ddply(y5,c("Region","interact"),transform,v2=quantile(Value,0.02))
y7 <- ddply(y6,c("Region","interact"),transform,v5=quantile(Value,0.05))
y8 <- ddply(y7,c("Region","interact"),transform,v10=quantile(Value,0.1))
y9 <- ddply(y8,c("Region","interact"),transform,v20=quantile(Value,0.2))
y10 <- ddply(y9,c("Region","interact"),transform,v33=quantile(Value,0.3333))
y11 <- ddply(y10,c("Region","interact"),transform,v66=quantile(Value,0.6666))
y12 <- ddply(y11,c("Region","interact"),transform,v80=quantile(Value,0.8))
y13 <- ddply(y12,c("Region","interact"),transform,v90=quantile(Value,0.9))
y14 <- ddply(y13,c("Region","interact"),transform,v95=quantile(Value,0.95))
y15 <- ddply(y14,c("Region","interact"),transform,v98=quantile(Value,0.98))
y16 <- ddply(y15,c("Region","interact"),transform,v99=quantile(Value,0.99))
z <- regionlabel(y16)

z_export <- y16[y16$GCM=="ge2m" & y16$SYR==yearname & y16$N=="N1" ,c(2,4,5,10,11,12,13,14,15,16,17,18,19,20,21)]

filenamecsv <- paste('../output/csv/',varname,'_',yearname,'.csv', sep = '')
write.table(z_export, filenamecsv, quote=F, col.names=F, append=T,sep=",")

#g <- ggplot() + 
#  geom_bar(data=z_export, aes(x=Region,y=v1)) +
#  ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
#  ylab(labelname) + xlab("") +	MyThemeLine

#plot(g)

#ggsave(g1,filename=name, dpi=250,width=10, height=6)

}
#----End --- PPOP_UNSH & PCAL: Bars across regions ----




